Novel Natural Inhibitors Targeting Enhancer of Zeste Homolog 2: A Comprehensive Structural Biology Research

The enhancer of zeste homolog 2 (EZH2) is a methylated modification enzyme of Histone H3-Lys 27. The high expression of EZH2 in cells is closely related to the progression, invasion, and metastasis of neoplasm. Therefore, this target is gradually becoming one of the research hot spots of tumor pathogenesis, and the inhibitors of the EZH2 enzyme are expected to become new antitumor drugs. This study used a series of virtual screening technologies to calculate the affinity between the compounds obtained from the ZINC15 database and the target protein EZH2, the stability of the ligand–receptor complex. This experiment also predicted the toxicity and absorption, distribution, metabolism, and excretion (ADME) properties of the candidate drugs in order to obtain compounds with excellent pharmacological properties. Finally, the ligand–receptor complex under in vivo situation was estimated by molecular dynamics simulation to observe whether the complex could exist steadily in the body. The experimental results showed that the two natural compounds ZINC000004217536 and ZINC000003938642 could bind tightly to EZH2, and the ligand–receptor complex could exist stably in vivo. Moreover, these two compounds were calculated to be nontoxic. They also had a high degree of intestinal absorption and high bioavailability. In vitro experiments confirmed that drug ZINC000003938642 could inhibit the proliferation and migration of osteosarcoma, which could serve as potential lead compounds. Therefore, the discovery of these two natural products had broad prospects in the development of EZH2 inhibitors, providing new clues for the treatment or adjuvant treatment of tumors.


INTRODUCTION
EZH2, namely, enhancer of zeste homolog 2, is a pivotal member of epigenetic regulatory factor Polycomb group (PcG) proteins. PcG proteins can lead to gene suppression through methylation modification (1), which comprises several essential molecules like Polycomb repressive complexes (PRCs). PRCs have inherent histone methyltransferase (HMTase) activity, which can inhibit gene expression through core histone methylation (2). PRC2 is of vital importance in PcG proteins, as it plays a role in the development of cancer (3). PRC2 consists of three subunits: EZH2, SUZ12, and EED, of which EZH2 and chaperone proteins are essential to correctly coordinate differentiation and proliferation of cells (4).
EZH2 has methyltransferase activity and can catalyze the methylation of histone H3-Lys 27 (H3-K27); it is essential for PRC-mediated gene suppression (5). Research had reported that human EZH2 was upregulated in different kinds of tumors like breast cancer, prostate cancer, and osteosarcoma (OS) (6). Cyclin-dependent kinase 1 (CDK1) promotes EZH2 ubiquitination by mediating the phosphorylation of Thr-345 and Thr-487 (T345 and T487) sites of EZH2 (7). And the posttranslational modifications of EZH2 are essential to improve its protein stability that related to the function of tumor cells and tumor metastasis, which could further lead to the accumulation of EZH2 and the occurrence of cancers (8,9).
In summary, EZH2 is related to different kinds of neoplasms, which was abnormally expressed and could serve as a therapeutic target (10)(11)(12). Therefore, inhibition of EZH2 protein could provide new ideas and methods in the treatment of cancers. GSK126 is a new type of competitive inhibitor targeting EZH2, which had begun tests in clinical trials (13). GSK126 significantly reduces the level of H3K27me3 in tumor cells by inhibiting the methyltransferase activity of EZH2, thereby inhibiting the growth of tumor cells such as human tongue squamous cell carcinoma and multiple myeloma cells (14,15). In addition, the appropriate concentration of GSK126 could also induce tumor cell apoptosis through the mitochondrial pathway (16). Research also reported that EZH2 may promote tumor invasion and metastasis by downregulating downstream targets such as Ecadherin and vascular endothelial growth factor (VEGF)-A (17,18). VEGF-A is an important cytokine that regulates angiogenesis, which is closely related to tumor metastasis (19). Research on EZH2 inhibitors has become hot spots in recent years, which has changed the treatment scheme as well as ideals dramatically. Nevertheless, novel efficient inhibitors targeting EZH2 still remained less. Consequently, more inhibitors regarding EZH2 were needed to discover from a natural medicine library in order to screen novel natural lead compounds and provide new clues in the discovery of EZH2 inhibitors. Existing studies had confirmed that EZH2 was highly expressed in OS patients and could serve as potential biomarker (11), while research on targeted therapy of OS targeting EZH2 had hardly been reported. Up to now, the research on EZH2 inhibitor GSK126 had made notable progress in different kinds of cancers, including prostate cancer cells and gastric cancer cells (20,21). Consequently, this study chose GSK126 as the reference compound to compare the pharmacological properties of the candidate compounds in order to discover more potential lead compounds targeting EZH2. Besides, this study aimed to validate whether EZH2 could serve as a therapeutic target in the treatment of OS.
Recently, natural products and natural extracts may be highly available compounds with proper biological activity that has potential medicinal value. They are therefore important sources for discovering, designing, and improving new drug skeletons (22,23). Extensive investigations have shown that natural products and their derivatives are currently playing an important role in the medical industry. It is already determined that the natural compounds from natural product database have considerable pharmacological potential (24,25). The second part of this study provides a theoretical basis and guidelines for discovering new inhibitors from natural product repository by screening inhibitory compounds related to EZH2. Besides, absorption, distribution, metabolism, and excretion (ADME) and toxicity prediction, ligand binding research, and molecular dynamics (MD) simulation were carried out on the selected candidate compounds, laying the foundation for the improvement of tumor drugs.

Discovery Studio Software and Ligand Library
The LibDock module of Discovery Studio 4.5 software (BIOVIA, San Diego, CA, USA) is used to screen better energy-optimized natural products, and the ADMET module can be applied to ADME analysis and the prediction of carcinogenicity and Ames mutagenicity. The CDOCKER module can be used to analyze the binding force between the products and the corresponding target of the protein and analyze the stability of the complex. This experiment selects the natural product database in the ZINC database to screen EZH2 inhibitors. The ZINC15 database is provided by the Irwin and Shoichet Laboratories in the Department of Pharmaceutical Chemistry at the University of California, San Francisco (UCSF), which provides a free virtual screening database of commercially available compounds.

Structure-Based Virtual Screening Using LibDock
The LibDock module was widely used in the drug development process (26). The LibDock module used a grid placed in the binding site and used polar and non-polar probes to calculate protein hot spots, then further used hot spots to arrange the ligands to form favorable interactions. Moreover, the study also used the Smart Minimiser algorithm and the CHARMm force field (Cambridge, MA, USA) to minimize the ligands (27). Then, all ligands' positions were adjusted and ranked according to the calculated ligand scores. The 2.5-Å crystal structure of human EZH2 [Protein Data Bank (PDB) identifier: 5WF7] and the structure of the inhibitor GSK126 were downloaded respectively from the PDB and ZINC15 database, and they were imported into the working environment of LibDock. Figure 5 shows the chemical structure of EZH2. The protein was prepared through several steps, including removing the crystal water and other surrounding heteroatoms, then adding hydrogen, protonation, ionization, and energy minimization. Among them, the energy minimization was realized by the CHARMM force field and the Smart Minimiser algorithm. In the case that the root mean square gradient tolerance was 0.1, the minimization performed 2,000 steps. After calculation, the binding site of the prepared protein was defined through the "Edit binding site" option. Analyzing the binding site of the ligand (GSK126) to generate the active binding site for docking the ligand with the receptor. Virtual screening was performed by docking the ligand exported from the database with the defined active binding site through the LibDock module. All compounds were grouped and ranked according to their LibDock score.

Prediction of Absorption, Distribution, Metabolism, and Excretion and Toxicity
The ADMET module of DS4.5 was used to estimate the adsorption, distribution, metabolism, and excretion properties of compounds. And the TOPKAT module of DS4.5 was employed to predict the carcinogenicity, Ames mutagenicity, and developmental toxicity potential in rodents. These pharmacological properties are fully considered when screening suitable EZH2 inhibitors to ensure the safety of the drug.

Molecule Docking and the Prediction of Drug Affinity
The CDOCKER module of DS4.5 was used for molecular docking research. CDOCKER is a tool to calculate high-precision docking results based on the CHARMM force field. During the docking process, the structures of the ligands are allowed to bend, while the structure of the receptor remains rigid. The CHARMM energy and interaction energy of each posture generated are calculated to reflect the binding affinity of the ligand and the receptor. Since crystal water molecules may affect the formation of receptorligand complexes, fixed water molecules are usually removed during the semi-flexible and rigid docking process, and hydrogen atoms are added to the protein to ensure the accuracy of the experiment (28,29). The crystal structure of EZH2 was obtained from the protein database, and the three-dimensional structures of ZINC3938642 and ZINC4217536 were obtained from the ZINC15 database. In order to verify the reliability of the results, this experiment also downloaded the reference compound GSK126 from the ZINC15 database. Similarly, the GSK126 was docked with EZH2 to calculate the root mean square deviation (RMSD) of the molecular docking conformation and compared it with the RMSD of the conformations of the ligandreceptor complex that are selected in this experiment. The binding site of EZH2 is defined as an area within a 5-Å radius from the geometric center of the ligand GSK126. In this experiment, the selected ligand was allowed to bind to the protein group residues in the binding site sphere. The identified hit structures were prepared and docked with the binding site of EZH2. Based on the numerical values of CDOCKER interaction energy, the different postures of each ligand-EZH2 receptor complex were generated and analyzed in detail.

Molecular Dynamics Simulation
Among the various postures predicted by the molecular docking program, the best binding conformation of the EZH2-inhibitor complex is selected as the object for MD simulation. The ligandreceptor complex is placed in an orthogonal box and solved with an explicit periodic boundary solvated water model. At the same time, to simulate the physiological environment, sodium chloride with an ionic strength of 0.145 was added to the system. Then, the system is subjected to the CHARMM force field and is relaxed through energy minimization (1,000 steps of steepest descent and 1,000 steps of the conjugated gradient). The reaction system was slowly driven from the initial temperature of 50K to the target temperature of 300K, the driving time was 2 ns, and the equilibrium simulation was performed when the time was 1 ns. The time for MD simulation (production) is 40 ns, and the time step is 2 fs. The simulation adopts the NPT (normal pressure and normal temperature) system at a constant temperature close to 300K, and the results were stored at a frequency of 0.02 ns. The Particle Mesh Ewald algorithm was used to calculate the long-range static electricity, and the linear constraint solver algorithm was used to fix all bonds involving hydrogen. Set the initial complex as the reference object. Use the Discovery Studio 4.5 analyze trajectory protocol to determine the structural properties, RMSD, and potential energy of the trajectory simulated by MD.

Cell Lines and Reagents
Human OS cell lines MG-63 (CL-0157), HOS (CL-0360), and human normal liver cell line LO2 (CL-0111) were purchased from Procell Life Science & Technology Co., Ltd. These cell lines were incubated in high-glucose Dulbecco's modified Eagle's medium (DMEM; Procell, Cat. no. PM150210), containing 10% fetal bovine serum (FBS; Gemini, USA) and 100 units/ml penicillin and 100 mg/ml streptomycin under normal cell culture conditions (37°C and 5% CO 2 ). Drug ZINC000003938642 was provided by Selleck Chemical Co. (Cat. no. S3668). The drug was dissolved in ultrapure water based on manufacturer-provided instructions to obtain the stock solution. Dimethylsulfoxide (DMSO) was not used to dissolve the drug in this study so that the toxicity effect on cells was negligible. Then, appropriate culture medium was added into the stock solution to configure different concentrations of the drug.

Cell Counting Kit-8 Assay
The standard Cell Counting Kit-8 (CCK-8) assay (provided by ApexBio, USA) was conducted to measure the cellular viability and proliferation of OS cells (HOS and MG-63) and human liver cell (LO2). Cell lines were plated into 96-well culture plates with a density of 3,000 cells/well overnight. Cells were treated with different doses of drug ZINC000003938642 for 24 h. The concentration gradients of each treatment were 0, 5, 10, 20, 40, 80, 160, 320, and 540 mmol/L. Cells were cultured for 1.5 h after addition of 10 ml/well CCK-8, and then the OD value of each well was measured based on the wavelength of 450 nm according to the microplate reader (BioTek Instrument, Synergy H1, USA).

Colony Formation Assay
Colony formation assay (CFA) assay was performed to detect the effects of different doses of drug on proliferation of tumor cells.
HOS and MG-63 cells were incubated into six-well plate with the density of 600 cells/well. After 24 h in culture, we configured cell culture medium with drug concentration of 100, 250 mmol/L; DMSO was not used in this study so the influence of DMSO on cells could be neglected. After 10 days of cultivation, the developed colonies were rinsed with phosphate buffered saline (PBS) twice and fixed in 4% paraformaldehyde, then 0.5% crystal violet solution was used to stain the colonies for half an hour. Lastly, we counted and described colonies according to microscopic examination.

In Vitro Scratch Assay
OS cells (HOS and MG-63) were cultured in six-well plate to assess the effects of drug on the migration of tumor cells. When the degree of fusion reached 90%, a 1-ml pipette tip was used to make a consistent cell-free area. Then, PBS was used to rinse twice to wipe off the cell debris, and serum-free medium was changed to culture, and different concentrations (0, 25, 50, 100, 250 mmol/L) of drug were used to treat cells and observe the scratch width at 0, 6, 12, 24 h. After corresponding time, we captured images of scraped area with phase contrast microscopy and measured the wounds and scratch width. The migration rate of OS cells was calculated as:

Western Blotting
OS cell lines (HOS and MG-63) were seeded into T25 culture flask and treated with different doses of drug ZINC000003938642 for 48 h. Then, proteins were extracted by radioimmunoprecipitation assay (RIPA), and bicinchoninic acid (BCA) protein assay was conducted to define protein standard curve and detect the protein concentration of each sample. Ten percent sodium dodecyl sulfatepolyacrylamide gel electrophoresis (SDS-PAGE) was used to separate proteins of samples, and then proteins were transferred to polyvinylidene difluoride (PVDF) membranes. Five percent nonfat milk dissolved in Tris-buffered saline and tween 20 (TBST) buffer was used to block the membranes for 2 h, after that, the membranes were incubated with primary antibodies [EZH2, c-Myc from Abcam and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) from Proteintech] at 4°C overnight. On the second day, the membranes were washed with Tris-buffered saline and tween 20 (TBST) three times and then horseradish peroxidase-conjugated secondary antibody was added to incubate the membranes for 1 h at room temperature. The membranes were visualized with enhanced chemiluminescence reagents to detect corresponding protein signals. Viber Bio Imaging tools were used to measure the band densities.

Apoptosis Assay
OS cells (HOS and MG-63) in log growth phase were inoculated into six-well plate and were treated with different concentrations of the drugs. After culturing for 24 h, cells were extracted through trypsin (without EDTA) and Annexin-fluorescein isothiocyanate (FITC)/propidium iodide (PI) double staining was performed according to the manufacturer's instructions. Lastly, the stained cells were analyzed by flow cytometry techniques; the apoptosis rates were examined by ACEA NovoCyte flow cytometry.

Pharmacophore Predictions of the Ideal Lead Compounds
After initial validation of the antitumor effects of the selected compounds, this study further analyzed their pharmacophore characteristics. Pharmacophore predictions of compounds were performed according to 3D-QSAR pharmacophore generation module, which generated up to 255 fits per molecule to represent a small molecule, and only fits with energy values within the energy threshold of 10 kcal/mol were finally preserved.

EZH2 Expression in Third-Party Database
To figure out the expression situation of EZH2 in OS, this study analyzed the expression values of EZH2 between normal and OS patients in Gene Expression Omnibus (GEO, https://www.ncbi. nlm.nih.gov/geo/) database. In total, three GSE series were analyzed including GSE14359, GSE33382, and GSE126209. As shown in Figure 1, results demonstrated that the expression of EZH2 in OS patients was significantly upregulated compared with that in normal patients (P < 0.05, Wilcoxon nonparametric test).

Fast Virtual Screening of Potential Inhibitors of EZH2
The SAL/SET domain of EZH2 protein is regarded as an important regulatory site for its enzymatic activity. Inhibitors bind to the SAL/SET domain of EZH2 by inserting into the ligand pocket of EZH2 and exerts the function of inhibiting the activity of EZH2: The small molecules binding to this site can prevent S-adenosyl methionine (SAM) from providing EZH2 with the methyl group needed to methylate H3K27me3, thereby reducing the enzymatic activity of EZH2. After SAM loses its methyl group, it is metabolized and hydrolyzed into intermediate products including S-adenosyl-L-homocysteine (SAH) and adenosine. S-adenosylmethionine is a methyl donor for onecarbon unit metabolism in organisms, and by moderately promoting the metabolic level of SAM, the activity of EZH2 can be inhibited (30). Based on this mechanism, inhibitors of enzyme activity against EZH2 could be identified. Therefore, this domain was chosen as the docking site for screening. The crystal structure of EZH2 was displayed in Figure 2, which contained the binding site sphere for docking, as well as the Ramachandran diagram of the protein, to check the rationality of EZH2 structure. Firstly, LibDock module of DS4.5 was performed to virtually screen small molecules that functioned in binding with the receptor protein EZH2. Downloading commercially available natural compounds from the ZINC15 database, a total of 13,537 ligands were generated by virtual screening. At the same time, the effective selective inhibitor GSK126, which could inhibit the activity of EZH2, was selected as the reference compound. After screening, 669 compounds were found with higher LibDock scores than GSK126 (LibDock score: 132.143). The top 20 compounds were listed in Table 1.

Absorption, Distribution, Metabolism, and Excretion Characteristics and Toxicity Prediction
By using the ADME and TOPKAT prediction module, we obtained the candidate 20 kinds of ligands and GSK126's pharmacological properties, including penetration of the blood-brain barrier, degree of human intestinal absorption, water solubility level, inhibitory effect on cytochrome P450 2D6, hepatotoxicity, and plasma protein binding properties ( Table 2). The water solubility prediction showed that 18 compounds could be dissolved in water relatively well. Among them, 10 compounds had a high solubility level (scores >2), which were greater than the reference compound GSK126 (moderate solubility, score: 1). For the degree of human intestinal absorption, 19 compounds had a good absorption effect, the same as GSK126, and ZINC000085826837 had a medium absorption level. Besides, seven compounds and GSK126 could be bound strongly by plasma proteins, while the remaining compounds did not have tight binding affinity and strong interactions with plasma proteins. Cytochrome P4502D6 (CYP2D6) was a key enzyme in the process of drug metabolism. Compounds involved in the screening had no inhibitory effect on CYP2D6. GSK126 was also predicted to be a non-inhibitor of CYP2D6. For liver toxicity, 12 compounds were predicted to be nontoxic drugs, while the remaining compounds and GSK126 were toxic to the liver. Subsequently, this experiment also calculated the safety properties of the candidate compounds and GSK126 through the TOPKAT module, including Ames (Ames mutagenicity), developmental toxicity potential (DTP), and rodent carcinogenicity [based on the United States National Toxicology Program (NTP) data set]. Experimental results displayed that those 12 compounds were non-mutagenic in long-term effect. It was predicted that four compounds were non-carcinogens and three compounds had no developmental toxicity potential. In addition, the reference compound GSK126 also predicted with pretty characteristics on Ames and NTP carcinogenicity, while it was computed with probability of DTP. The detailed information of the indicators among compounds and GSK126 were shown in Table 3. Based on the above data, ZINC000004217536 and ZINC00000 3938642 were neither CYP2D6 inhibitors nor hepatotoxicity drugs. Moreover, they were predicted to be free of Ames mutagenicity and rodent carcinogenicity. Consequently, ZINC000004217536 and ZINC000003938642 were analyzed to be candidate drugs with  high safety and were selected for further study. The detailed chemical and crystal structures of these compounds were shown in Figure 3.

Ligand-EZH2 Binding Analysis
In order to study the binding mechanisms between the ligand and receptor EZH2, CDOCKER module was conducted to dock ZINC000004217536, ZINC000003938642, and GSK126 at the regulatory site of EZH2, and the corresponding CDOCKER potential energy of these complexes was calculated, as shown in Table 4. Hydrogen bonds and p-p interactions between EZH2 and these compounds were analyzed (Figures 4, 5).  Table 5.

Molecular Dynamics Simulation
The best binding conformations of each compound-EZH2 complexes (ZINC000004217536-EZH2 and ZINC000003938642-EZH2) were obtained from precise docking program CDOCKER and applied for the following experiment. In this study, the stability of the ligand-EZH2 complex under in vivo circumstance was predicted by MD simulation module. The predicted results were shown in Figure 6, including energy values ( Figures 6A, D) and RMSD curve ( Figures 6B, E) diagram of each complex. The orbitals of all complexes reached equilibrium after 100ps. The complexes' RMSD and energy values like total energy, potential energy, and electrostatic energy all got stabilized over time. Results suggested that hydrogen bonds formed by the compound and EZH2 and the p-p-related interactions contributed a great effect on the stability of these complexes. Furthermore, chemical bonds heatmap after MD also illustrated that these chemical bonds, which contributed largely to the stability of complexes, could still exist steadily with the progression of MD in the body (Figures 6C, F). Based on the above evaluation indicators, the complexes formed by  ZINC000004217536 and ZINC000003938642 with EZH2 could exist stably in the internal environment. Consequently, these two compounds could interact with EZH2; they also had a regulatory effect on EZH2 like the reference compound GSK126 did.

ZINC000003938642 Reduced Proliferation of Osteosarcoma Cells
To test the antitumor effects of compounds screened in this study, we purchased one of the two compounds, ZINC000003938642, for    Figures 7A, B). Subsequently, to validate the toxicity of drug to liver cells, LO2 cell line was conducted and measured by CCK-8. Results revealed that drug ZINC00000 3938642 did not inhibit the proliferation of human normal liver cells in a dose-dependent manner and time-dependent manner, which still had a high cellular viability even when subjected to the highest dose ( Figure 7C). We then performed CFA to further determine the antitumor effects of drug in OS cells. As shown in Figure 7D, after 10 days of cultivation with different drug concentrations (100 and 250 mmol/L), both HOS and MG-63 cells showed fewer and smaller clonogenicities in Petri dishes with drug group than with the control group. The numbers of clone formation in drug groups were significantly lower than those in control groups (P < 0.05) (Figures 7E, F).

ZINC000003938642 Inhibited Migration of Osteosarcoma Cells
To analyze the effects of drug ZINC000003938642 on OS cell migration, scratch assay was further performed. The width of scratched areas was measured at 0, 6, 12, and 24 h of scratch, and the scratch width represented the migration capacity of OS cells (Figures 8A, D). As shown in Figures 8B, C, results exposed that with the extension of time and increase of drug concentration, the migration rate of OS cells to the scratch area decreased significantly (P < 0.05).

ZINC000003938642 Induced Apoptosis in Osteosarcoma Cells
Flow cytometry and Annexin-FITC/PI double staining were used for apoptosis assay to detect the effects of drugs on programmed cell death. The apoptotic rates of HOS and MG-63 cells were detected after being treated with concentrations of ZINC000003938642 (0, 50, 100, and 250 mM) for 24 h. As shown in Figure 9A, results illustrated that the apoptotic rates increased with the increase of drug concentration both in HOS and MG-63 cells. Consequently, live cells were predominant in the control (0 mM) groups, while apoptotic cells were predominant in the drug-treated groups ( Figure 9B; P < 0.05).

ZINC000003938642 Reduced the Expression of EZH2 and Its Downstream Gene C-Myc
The expression of EZH2 was detected by Western blot analysis. As shown in Figures 10A, C, the expression level of EZH2 was inhibited by ZINC000003938642 in both HOS and MG-63 cells, and its inhibitory effect displayed a dose-dependent manner (P < 0.05; Figures 10B, D). c-Myc was the downstream target gene of EZH2, and as an oncogene, the expression levels of c-Myc were also reduced by ZINC000003938642 in both HOS and MG-63 cells ( Figures 10A, C), and its inhibitory effects also demonstrated a dose-dependent manner ( Figures 10B, D).

Ligand Pharmacophore Predictions
After initially verifying the antitumor effects of the candidate compounds, this study further analyzed the pharmacophore characteristics of these two compounds in order to observe the potential modification site on compounds. As shown in Figures 11A, B, computational results illustrated that there were 51 features in ZINC000004217536 and 69 features in ZINC000003938642, among which, ZINC000004217536 possessed 18 hydrogen bond acceptors, 21 hydrogen bond donors, five hydrophobic centers, one ionizable positive, and one ring aromatic. As for ZINC000003938642, it possessed 22 hydrogen bond acceptors, 38 hydrogen bond donors, two hydrophobic centers, five ionizable positive, and two ring aromatics.

DISCUSSION
Statistics published in these years show that malignant tumors are still the main causes of death among residents in many countries   (31). OS is one of the most common malignant tumors of mesenchymal tissues, occurring mostly in the metaphysis of long bones in adolescents; it is a differentiation-defective disease caused by dysdifferentiation of osteoblasts and/or epigenetic changes (32). The harm of malignant tumors to humans is not only a threat to the lives of patients but also the physical pain, mental pressure, and the economic burden they bring to patients (33). EZH2, an essential component of the epigenetic regulatory factor PcG and a catalytic subunit of PRC2, is involved in regulating the methylation of lysine 27 (H3K27) of histone H3 and is highly expressed in a variety of tumors. It plays an important role in regulating gene transcription and gene silencing and participates in the growth, proliferation, and metastasis of tumor cells.
In recent years, EZH2 has become a popular target for cancer therapeutics, and the research of EZH2 inhibitors and their combined application with other antitumor drugs in clinical practice has broad prospects (34). However, relatively few inhibitors of EZH2 had been discovered and analyzed. Existing research had reported the high expression of EZH2 in OS patients. Currently, few studies have been conducted on the efficacy of EZH2 inhibitors in OS. GSK126, which is relatively a mature EZH2 inhibitor, was applied in this study to analyze the  antitumor effect and molecular mechanism on OS, and it was regarded as the reference compound to compare pharmacologic properties with novel ligands. Although GSK126 has certain antitumor functions, it still has some limitations. Relevant studies have shown that EZH2 can produce drug resistance through allelic mutations and protein conformation changes (35,36). For the purpose of overcoming the drug resistance of EZH2, it is necessary to develop new inhibitors. Furthermore, GSK126 had low solubility and rapid plasma clearance, resulting in low bioavailability, and GSK126 had hepatotoxicity, which leads to unsatisfactory effects of highdose GSK126 in the process of inhibiting EZH2.
In this study, we aimed to discover more potential lead compounds of EZH2. The available natural compound structures were downloaded from the ZINC15 database for virtual screening, then ADME, TOPKAT, CDOCKER, and other modules were applied to perform ADME prediction, rodent carcinogenicity and Ames mutation prediction, ligandreceptor binding studies, and MD simulations. The LibDock scores suggested the degree of energy optimization and conformational stability between compound and receptor. Compounds with higher LibDock scores illustrated better energy optimization and more stable conformation than compounds with lower scores. Calculation results of the LibDock module showed that a total of 13,537 compounds could be stably combined with EZH2 after fast docking method. Among these ligands, 669 compounds had higher LibDock scores than the reference compound GSK126 (LibDock score: 132.143), suggesting that the stability and energy optimization effect of the complex substances formed by these 669 compounds with EZH2 were more stable than GSK126-EZH2 complex. Based on ranking of the LibDock scores, the top 20 compounds with the highest scores were screened out and tested in next steps. ADME and toxicity prediction were performed to evaluate the pharmacological properties of these selected compounds. After analysis, results elucidated that ZINC000004217536 and ZINC000003938642 had satisfactory intestinal absorption capacity, and these two compounds had no obvious inhibitory effect on CYP2D6, no hepatotoxicity, and low binding affinity property with plasma protein, which suggested the good selectivity of these drugs; they could avoid rapid clearance by plasma so as to behave the best pesticide effect. Furthermore, compared with other compounds, ZINC000004217536 and ZINC000003938642 were predicted not to have Ames mutagenicity and rodent carcinogenicity, and they had less developmental toxicity potential. Consequently, they were considered as ideal candidate compounds with pharmacologic properties and higher security in the body; these characteristics were enough to be considered as the most potential lead compounds. Based on the above results, ZINC000004217536 and ZINC000003938642 were reasonably recognized as high-quality medicinal materials; these two compounds had broad application prospects in drug development and design. Although other drugs on the list had certain negative effects such as developmental toxicity and Ames mutagenicity, other pharmacological properties were relatively moderate, so they also had a certain potential in drug improvement, which could be achieved by adding or deleting specific functional groups or atoms to reduce their negative effects. In summary, it was determined that ZINC000004217536 and ZINC000003938642 were the most potential lead compounds, and more analyses were further performed. Subsequently, we analyzed chemical bonds and the binding mechanisms between candidate compounds, GSK126 and EZH2. Precise docking method CDOCKER module was conducted; results showed that the CDOCKER interaction energy of ZINC000004217536 and ZINC000003938642 was significantly lower than that of the reference ligand GSK126, proving that the affinity of these two compounds with EZH2 was higher than GSK126-EZH2 in real situations. After that, through molecular detection analysis of their chemical structures among these complexes, results illustrated that the chemical bond and interaction force of the complex formed by EZH2 and the two candidate compounds were stronger, which further explained that they may have a competitive inhibitory effect on the regulatory site of EZH2 and thus inhibit the activity of EZH2, finally producing antitumor effects.
Ultimately, MD simulations were performed to predict the stability of the complexes formed by the candidate compounds and EZH2 in the internal environment. By calculating the RMSD and energy values of these ligand-EZH2 complexes, the RMSD curve and energy curve were drawn. Results showed that the trajectories of the complexes reached equilibrium after 100ps, and the RMSD and energy values of these complexes tended to be stable over time, indicating that the two complexes could exist stably in the internal environment. Furthermore, chemical bonds heatmap elucidated that these chemical bonds, which contributed remarkably to the stability of the complex, could keep steady with the progression of the MD. Consequently, the compounds selected in this study bonded tightly to EZH2, and they were capable of existing stably in the body, thereby exerting corresponding pharmacological functions. Therefore, they have great potential in the development of EZH2 inhibitors. It is noteworthy that the reference compound GSK126 chosen in this study served as a known effective synthetic EZH2 inhibitor; the effects of natural compounds were hardly better than GSK126 in vitro or in vivo. The role GSK126 played in this process was to provide a primitive crystal complex for us to compare the binding mode and give us an initial active binding sphere, and it did not provide any guiding significance for the comparison between GSK126 and ZINC000003938642 in antitumor aspects.
Currently, existing studies pointed out EZH2 could serve as a therapeutic target regarding OS (11), while few studies focused on targeted therapy of OS targeting EZH2. Consequently, this study preliminarily discussed the effects of newly found compounds against OS. To prove the pesticide effects of our newly found compounds against OS and the reliability of the screening method in this study, we selected one of the candidate compounds, ZINC000003938642, and performed a series of in vitro experiments including CCK-8, CFA, scratch assay, Western blot, and apoptosis assay. In CCK-8 assay, results pointed that the cellular viability in OS cells had a dosedependent decrease when treated with drug ZINC000003938642, while the drug was relatively well tolerated for human liver cells LO2. This finding implied that this drug was relatively nontoxic in term of hepatotoxicity, which was also consistent with our predictions in structural biology part that ZINC000003938642 was a nontoxic drug. In CFA, the numbers and size of clonogenicities in drug group were significantly less than those in control group in both HOS and MG-63 cell lines, which was consistent with results that the proliferation of OS cells was reduced by drug in CCK-8 assay and that the effects were dosedependent. Scratch assay revealed that the wound area in control group decreased more sharply than that in drug group with time. As for apoptosis assay, flow cytometry results visualized that the percentage of apoptotic cells increased with the drug increasing, the apoptotic rates of HOS and MG-63 cells treated with high drug dose groups were significantly higher than those of control (0 mM) group (P < 0.05). Western blot analysis revealed that EZH2 expression decreased with increasing drug concentrations. Since c-Myc is the downstream target of EZH2, and as an oncogene (37)(38)(39), the expression level of c-Myc could also reflect the inhibitory effects of the drug on EZH2. Results displayed that the downstream oncogene c-Myc was also inhibited by the drug in a dose-dependent manner, implying that drug ZINC000003938642 could serve as a potential EZH2 inhibitor. These experiments suggested the ability of drug to inhibit the proliferation, migration, and EZH2 and c-Myc expression of OS cells, which indicated that drug ZINC00000 3938642 found in this study was an effective inhibitor regarding OS, and EZH2 was a therapeutic target against OS. The screening of ideal lead compounds is a key step in drug design and development. Regarding the pharmacophore predictions of ZINC000004217536 and ZINC000003938642, they possessed a number of pharmacophores, which elucidated that based on these skeletons of these two compounds, the modification and refinement of the drug could be conducted to further make a whole new design. The natural compounds discovered in this study are of great significance in the development of EZH2 inhibitors. This study provided evidence for the targeted treatment of OS regarding EZH2 and may have the potential to provide better methods for tumor treatment. Besides, in the field of pharmacology, more research could be studied like modifying the molecular structure of the drugs to reduce the toxicity and mutation to continuously improve the pharmacological effect of the inhibitor.

CONCLUSIONS
This study used a series of virtual screening techniques and discovered two natural compounds, ZINC000004217536 and ZINC000003938642, which have the function of inhibiting the active subunit EZH2 of PRC2. These two compounds bind tightly to the target protein. Additionally, they have no carcinogenicity and toxicity, so they can be regarded as potential EZH2 inhibitors. In vitro experiments confirmed that drug ZINC000003938642 could inhibit the proliferation and migration of OS, which could serve as potential lead compounds. This study not only provided the pharmacological properties of candidate drugs but also provided meaningful materials for further research of EZH2-targeted inhibitors.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material. Further inquiries can be directed to the corresponding authors.